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Theoretical analysis and fully atomistic molecular dynamics simulations reveal a Brownian ratchet 
mechanism by which thermal fluctuations drive the net displacement of immiscible liquids confined 
in channels or pores with micro- or nanoscale dimensions. The thermally-driven displacement is 
induced by surface nanostructures with directional asymmetry and can occur against the direction 
of action of wetting or capillary forces. Mean displacement rates in molecular dynamics simulations 
are predicted via analytical solution of a Smoluchowski diffusion equation for the position probability 
density. The proposed physical mechanisms and derived analytical expressions can be applied to 
engineer surface nanostructures for controlling the dynamics of diverse wetting processes such as 
capillary filling, wicking, and imbibition in micro- or nanoscale systems. 


INTRODUCTION: NANOSCALE WETTING AND BROWNIAN RATCHETS 


Advances in nanofabrication and characterization techniques have enabled the engineering of nanostructured sur¬ 
faces with geometric features as small as a few nanometers w- At nanoscales, the interplay between intermolecular 
forces, Brownian motion, and surface structure can give rise to complex interfacial phenomena that are challenging 
for the application of conventional, continuum-based and deterministic, models mu. For example, nanoscale surface 
structures can induce energy barriers that lead to wetting processes governed by thermally-activated transitions be¬ 
tween metastable states (SHU]. These thermally-activated transitions can result in directed transport of fluids and 
solutes when there is directional asymmetry of the energy barriers induced by the physicochemical structure of the 
confining surfaces HMg. Analogous mechanisms for rectification of thermal motion into directed transport underlie 
fundamental biological processes such as selective charge transport in ion channels or translocation of proteins across 
cellular membranes. Physical systems where thermal fluctuations are able to drive net directional motion, while per¬ 
forming work against “load” or resistance forces, are known as thermal ratchets or Brownian motors and have been 
extensively studied in the framework of statistical physics mm- 

Thermal ratchets can operate without thermal or chemical gradients provided that the system has not reached all 
necessary conditions for thermodynamic equilibrium [201 [21]. A variety of novel nano/microfluidic devices perform as 
thermal ratchets to accomplish the handling, separation, and detection of diverse solutes (e.g., DNA, macromolecules, 
ionic species) and/or colloidal particles with an unprecedented precision [22l - (25] . These devices usually work with 
single-phase fluid solvents and must combine external electromagnetic fields, electrolyte solutes in proper concen¬ 
tration, and formation of electric double layers in order to induce energy landscapes with directional asymmetry 
(i.e., ratchet potentials). A different class of ratchet systems involving multiphase fluids has been demonstrated to 
produce “self-propulsion” of micro- or millimeter-sized droplets by combining micro/nanostructured surfaces, ther¬ 
mal/chemical gradients, and/or mechanical vibration pM3T] . Self-propulsion mechanisms in these multiphase systems 
are attributed to diverse dynamic phenomena, such as capillarity and contact angle hysteresis [28 , 29] , or evaporation 
flows and the Leidenfrost effect mm, where thermal fluctuations play a secondary role. 

There is a class of multiphase (two fluid) system that can perform as a thermal ratchet under isothermal and 
incompressible conditions, with or without the presence of electrolyte solutes and net surface charge. In this class of 
system the thermal ratchet mechanism is enabled by surface nanostructures that induce surface energy barriers with 
directional asymmetry. The particular configuration considered in this work, illustrated in Fig. [TJl, consists of two 
macroscopically immiscible liquids (fluid-1 and fluid-2) confined in a slit-shaped channel or pore of height h Q , length /, 
and width w <C h Q . The surfaces confining the fluids are chemically homogeneous and neutrally charged. One of the 
surfaces has a terraced structure with regular tread length s and riser height r [cf. Fig. FF of nanometric dimensions. 
Similar terraced structures have been synthesized on crystalline substrates via diverse nanofabrication techniques 
such as wet etching, high-temperature annealing, and deposition of epitaxial films HJEH33. The studied terraced 
structure with Ns steps reduces the local height of the channel according to h(x ) = h(0) — r[x/s\ for 0 < x < (Ns + l)s 
(here, [x\ = floor (x) is the floor function and x is the coordinate in the longitudinal direction). In the presence of an 
interface between two immiscible fluids, the interplay between thermal motion and surface energy barriers induced by 
the nanoscale structure can drive imbibition and filling/drainage processes in micro/nanoscale channels or pores for 
a range of wettability conditions unanticipated by conventional wetting models. 
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THEORETICAL DESCRIPTION: THERMALLY DRIVEN WETTING 


Analytical descriptions of thermally-driven wetting processes must consider that atoms or molecules in a liquid- 
fluid interface undergo thermal motion. We will analyze the case of unidirectional motion described by the average 
position x(t) of all atoms of the first fluid species (fluid-1) that he at the front liquid-liquid interface [cf. Fig. ER 
Adopting the average interface position to describe the dynamics of the confined molecular fluids implies projecting 
the (multidimensional) system energy landscape onto a one-dimensional profile U(x) along a “reaction coordinate” x. 
The sequence of random displacements of the front interface position can be statistically described by the conditional 
probability density p(x, t ) = p(x, t\x Q , t Q ); here, x Q = x(t Q ) is the average interface position observed at a time t Q . The 
stationary probability density p(x,t —)> oo) = Z -1 exp [—U(x)/k B T] is prescribed by the free energy profile U{x) and 
the thermal energy k B T] here, Z is the corresponding partition function, k B is the Boltzmann constant and T = const, 
is the system temperature. Assuming overdamped Brownian dynamics, the time evolution of the probability density 
p(x, t ) is governed by the Smoluchowski diffusion equation 


d_ 

dt 


p(x,t) 


d 1 dU~ 
dx 2 £(x) + dx £(x) dx 


P(x,t) 


( 1 ) 


where £(x) is the local friction coefficient or resistivity (i.e., the inverse of the mobility). For the studied conditions 
we consider a linear friction force Ff = —t^dx/dt that is mainly due to hydrodynamic effects and thus 


£(x) = k B ph(x)w 
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S 

(h Q — nr) 2 


( 2 ) 


where k B is a drag coefficient, fi is the shear viscosity of the confined fluids, and Ns is the total number terraces in the 
structure. For analytical simplicity, we consider in Eq. [2] the case that both fluid-1 and fluid-2 have the same viscosity 
/x = pi = /X 2 ; expressions for /xi 7 ^ /xi can be readily derived using similar hydrodynamic arguments. Analytical 
estimates of the drag coefficient k B in Eq. [2] can be obtained by making simplifying assumptions about the modes 
of translation of the interface and the hydrodynamic velocity profiles induced. A drag coefficient k B — 4 is obtained 
by naively assuming that a linear flow profile (i.e., Couette flow between two flat surfaces) develops after the sudden 
displacement of the contact line on one wall, while the contact line on the opposite wall remains stationary. 

For isothermal and incompressible conditions and assuming sharp interfaces, the free energy profile is determined 
by surface energy contributions 


U(x) — 7 [wh(x) — cosOy Ais(x)] + const. 


( 3 ) 


Here, 7 is the interfacial energy of the liquid-liquid interface, Oy is the Young contact angle measured on the fluid-1 
phase, and A\s is the area of the interface between the fluid-1 and solid phases. In the studied configuration the 
distance between the front and rear liquid-liquid interfaces is larger than the length of the terraced structure, which 
allows us to simplify the analysis. While the front interface moves inside each terrace there is a linear change in the 
interfacial area Ais(x) that results from conservation of volume, and thus we have 


dU(x) 

dx 


= —F n = — 2rry cos Oy 


wr 

h Q 


( 4 ) 


for ns < x < (n + 1 )s (n = 1 , Ns). 

The energy profile U(x) (Eq. [3| exhibits sharp energy increments AU- cs jwr when moving in the negative re¬ 
direction across the edge of a terrace at position x n = ns (n = 1, N$). Hence, periodic energy barriers AU- at 
regular steps 5 hinder backward random displacements as the liquid-liquid interface undergoes thermal motion along 
the terraced surface. The time to cross over an energy barrier induced by nanoscale surface features can be predicted 
via Kramers theory of thermally-activated transitions, as documented in prior work PJJEE3]. For the present analysis, 
it suffices to recognize that a mean time T_ oc exp (A U-/k B T) must elapse before observing a backward displacement 
of the interface over a terrace edge. Therefore, over a time interval t < T_ the presence of the energy barrier AU- can 
be treated as a reflective boundary condition for Eq. [l] imposed at the edge of each terrace. Furthermore, the friction 
coefficient (Eq. [2| and the wetting or capillary force (Eq. [4| remain constant within each of the Ns terraces. Hence, 
we have £ n = £(x) = const and F n = —dU(x)/dx = const for ns < x < {n + l)s, which facilitates the analytical 
solution of Eq. [l] In order to solve Eq. [l] within each terrace for which the edge lies at x n = ns , it is convenient 
to introduce the dimensionless position x = (x — x n )/s and time t = (t — t 0 )(k B T/^ n s 2 ). In addition, we introduce 
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the dimensionless force parameter C n = F n s/k B T given by the ratio of the work performed by the wetting force to 
the thermal energy within each terrace. For a reflective boundary at x n and initial condition x Q = x n the analytical 
solution of Eq. |T] is given by [35, 3d p n (x,t) = p(x,t; C n ) = p n /s, where 


PnW) 



(X + CrXf_ 

if 


+ ~c n exp [~C n x] erfc 


(x - C n t ) 

Vi! 


( 5 ) 


is a dimensionless function determined by the variables x > 0 and t > 0 for a given force parameter C n . The mean 
displacement of the front liquid-liquid interface according to Eq. [ 5 ] is (x(t) — x n ) = sX n {t ), where 

rOO 

Xn(t) = / p n (T, t)xdx (6) 

Jo 

is the dimensionless mean displacement within the n-th terrace. Taking the upper integration limit to infinity in Eq. [6] 
is an approximation valid for finite times t (n — N$ + 1 ) 2 s 2 k B T/^ n . From Eq. [6] we can estimate the mean time 


tn — tl T 


k B T 




k =1 


( 7 ) 


at which (x(t n )) = x n (n = 1 , TV's + 1). Here, X n 1 is the inverse of the dimensionless mean displacement function in 
Eq. [ 6 ] and t\ is the time at which (x(ti)} = s. 

A few comments aobout the derived expressions are in order. Analytical integration in Eq. [6] is feasible and thus 
the dimensionless time T n = (£ n +i — s 2 to traverse the n-th terrace (n = 1, Ns) can be obtained by solving 

the implicit equation X n {T n ) = E The inverse function Ae 1 (1) = 77 = T(C n ) in Eq. [ 7 ] is uniquely determined by 
the dimensionless force parameter C n = F n s/k B T ; the function T n 00 diverges at C n = — 1 and decays 7^ —>■ 0 for 
C n —>• 00 . An approximate explicit expression 77 = (4/7 tC^)( 1 — + (7r/2 )C n ) 2 can be obtained from a first-order 

Taylor expansion of Eq. [ 5 ] about C n = 0; this approximation yields less than 10% error for \C n \ < 1/2. For a neutral 
contact angle 0y = 90° the wetting force vanishes F n = 0 and thus 77 = 7r/4forn = l, N$. Moreover, under neutral 
wetting conditions Eq. [5] predicts a mean displacement (x(t)) = \j2D 0 t characteristic of diffusive processes, with an 
effective diffusivity D 0 = (2/7r)(k B T/^ n ). Notably, forward liquid displacements (x(t)) > 0 against wetting forces 
F n < 0 are expected to occur for contact angles 0y > 90° provided that — 271^ cos Oyw x ( r/h Q ) < k B T/s. 


MOLECULAR DYNAMICS SIMULATIONS 


In order to verify analytical predictions and underlying assumptions we perform fully atomistic molecular dynamics 
(MD) to simulate the dynamics arising from couplings between Brownian motion, hydrodynamic effects, and wetting 
forces. The MD techniques employed in this work are extensively described in the literature (33®! and prior work 
by the authors naiisi. The simulated system [see Fig. la] comprises two monatomic liquids (atomic species i = 1,2) 
labeled as fluid -1 and fluid- 2 , and a crystalline solid (atomic species i = 3). Atomic interactions are modeled by 
pairwise Lennard-Jones potentials ulj(v) = s[(r/a )~ 12 — Qj(r/cr) -6 ], where e is a characteristic interaction energy, 
Cij = Cji is a symmetric attraction coefficient between species (i,j = 1,3), r is the interatomic distance between any 
two atoms, and a is the diameter of the repulsive core, which roughly correspond to the atomic diameter. Following 
conventional techniques for computational efficiency, atomic interactions are not computed for r > 2.5 a. The time 
scale of atomic displacements is r = a^/m/e and is in the order of picoseconds for simple molecular liquids; the atomic 
mass m of the liquid species is set to be equal for both modeled liquids. The equations of motion are integrated using a 
fifth-order prediction-correction algorithm with a time step 5 = 0.004r. A Nose-Hoover thermostat 133 ESI HU models 
the interaction with a thermal bath, regulating the system temperature to a prescribed value. In all MD simulations 
in this work the prescribed temperature is T = 2s/k B T, the mean number density is (p) = 0.8/cr 3 and the shear 
viscosity is p = 2.4m/(crr) for both liquids. Solid atoms form a face-cubic-centered (fee) lattice with uniform spacing 
Ax = 0 . 8 - 1 / 3 cr; typical values of Ax for crystalline solids range between 0.2 and 0.5 nm. The set of values employed 
for the attraction coefficients (C 12 = 0.5, C 13 = 0.8, 0.75 < C 23 < 0.95, and c B = 1) render two macroscopically 
immiscible fluids with a liquid-liquid interfacial tension 7 = 1.2 e/a 2 and Young contact angles near neutral wetting 
conditions 85° < 0y < 100°. 
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(a) (b) 0Y = 90° s = 6Ax r = 2bx 



FIG. 1. Modeled system and geometric configuration, (a) Slit channel (height h Q , length /, and width w) confining two 
immiscible liquids. A terraced structure (tread s and riser r) of nanoscale dimensions lies on the bottom wall, (b) Volume 
fraction <p(t) in different MD realizations and its ensemble average (4>(t)) (time reported in MD units r ~ 1 ps). MD simulations 
correspond to s = 6 Ax and r = 6Ax under neutral wetting conditions Oy — 90°. Dashed lines indicate filling fractions 
(j) n ss (f)(x n ) for x{t) = ns (n — 1 ,Ns + 1). (c) Sequence of three time instances in MD simulations for contact angles Oy — 85- 
95°. Imbibition of fluid-1 into the volume V above the terraced structure is observed for 6y > 90° when capillary forces are 
neutral or negative F = —dU/dx < 0. 


RESULTS: THEORY AND MOLECULAR DYNAMICS 


The studied geometric configuration [cf. Fig. consists of a slit nanoscale channel of height h Q = 30-60Ax, length 
l = 70Ax, and width w = lOAx; periodic boundary conditions are applied in the x- and ^-direction. The bottom 
wall has different surface structures with Ns = 5-11 terraces of length s = 3-9Ax and riser height r = 1-2Ax. 
The volume confined above the terraced structure is V = ws Y^k=il^o ~ At]. The fluid-1 phase fills a fraction 0(x) = 
xw(h 0 -\-s)/V —x 2 (wr/2V s)-\-0{r/h(x)) of the volume V as the front liquid interface moves within s < x(t) < (Ns + l)s 
[cf. Fig. [lji]. In our MD simulations, the number density is constant and the filling fraction = n\/{n\ + 722 ) is 
directly computed from the number of atoms of fluid-1, n 1 , and fluid-2, 77 - 2 , occupying the volume V. At initialization, 
the volume V is fully occupied by the fluid-2 phase and thus </>(£ = 0) = 0. 

The filling fraction </>(£) for six MD realizations [42] and the ensemble average (0(t)) are reported in Fig. [TJd for 
one of the studied structures (s = 6Ax, r = 2 Ax) and neutral wetting conditions (cis = C 23 = 0 . 8 ). As seen in 
Figs.jTJ), the mean filling fraction (</>(£)) uniformly increases with time while individual MD realizations exhibit rapid 
transitions between “long-lived” metastable states at specific values <fi n ~ <j)(x n ) (n = l, Ns + 1)- The imbibition 
process occurs within a range of Young contact angles 0y ^ 90 [cf. FigJT^] for which capillary forces can be neutral 
or negative (F n = —dU/dx < 0). Displacement beyond the last terrace edge, where (j)(xjsr s + 1 ) = 1, is prevented by a 
large energy barrier AU+ = 'yNsrw. 

Theoretical predictions from Eq. [ 7 ] for mean displacements, (x(t n )} = x n , and filling fractions, (0(t n )) = (j)(x n ), 
are in close agreement with MD simulations reported in Figs. [2]-[3] when using kn = 4-5.5 in Eq. [2] for different 
structures. As seen in Fig. [ 2 ^ 1 , the mean filling rate is approximately constant (</>(£)) — const under neutral wetting 
conditions 0 = 90° for the different studied structures. Notably, mean displacements (x(t)) collapse to a unique curve 
when scaling by the terrace length s and diffusion time To = s 2 ^ 0 //c^T as showed in Fig. [2 ]d; here, £ 0 = £(0) is the 
resistivity in Eq. [ 2 ] for h(x = 0) = h Q . As shown in Fig. [ 2 ]:, increasing the channel height h Q for a given length l 
enhances the interface displacement rate by decreasing the local resistivity £(x). The effect of varying the contact 
angle 0y and thus the wetting force F n (Eq. [ 7 ]) can be seen in Fig. [2jl. For 0y < 90° the mean displacement is 
enhanced by positive wetting forces F n > 0. As predicted, positive displacements still can be observed for 0y > 90° 
when F n s/kBT > —1; the most adverse wetting force Fn s — —ksT/s corresponds to 0y cs 98° (s = s/Ax = 6 , 
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FIG. 2. Thermally-driven wetting on terraced nanostructures. Solid lines: linear interpolation between coordinate pairs (x n ,t n )] 
x n — ns (n = 1, Ns + 1) and t n given by Eq. [t| Markers: ensemble average quantities from MD simulations [see legends]. 
Geometric dimensions are in lattice units: channel height h = h 0 /Ax = 30-60, terrace length s = s/Ax = 3-6, riser height 
r — s/Ax = 1-2. (a)-(c) Imbibition under neutral wetting 0y = 90°: (a) mean filling fraction (</>(£)) (time in MD units 
r); (b) dimensionless displacement (x(t/Tn))/s where Td — Td — s 2 £o/&bT; (c) dimensionless displacement (x(t))/s for 
different channel heights, (d) Dimensionless displacement {x(t/Tu))/s for contact angles 0y = 85°, 90°, 95°, and 98° for which 
Cn s — 0-7, 0, -0.7, and -1, respectively. 



FIG. 3. Failure of the thermal ratchet mechanism, (a) Dimensionless displacement (xit/To))/s in MD simulations for large 
terrace length s — 9 and two riser heights s =1-2 (6y = 90°, h = 30). (b) Schematic of the failure and onset of diffusive motion 
with zero mean displacement for T_ < t < T_ + T n . 


r = r/ Ax = 2) in Fig. |2]i. 

Results reported in Fig. [2] confirm that a thermal ratchet mechanism mediated by surface energy barriers A U- — 
7 wr can lead to thermally-driven transport of immiscible liquids in slit channels or pores. This phenomenon is 
expected to occur provided that the time T_ oc exp (—^wr/ksT) to cross over each energy barrier is larger than the 
time t n - t n - 1 (Eq. 0 to traverse the n-th terrace. Indeed, MD simulations for large terrace length s = 9Ax and 
small riser height r = Ax [see Fig [3^i] report a decay in the mean displacement rate after a time T_ Increasing 

the energy barrier A U- by increasing the riser height to r = 2Ax prevents the observed decay in the displacement 
rate [see Fig [3^]. As illustrated in Fig [3p, the mean displacement rate is expected to vanish (x(t)} 0 for t > T_, 

after which unbiased Brownian motion persists over a time T n = s 2 £ >n /kBT = Td x (h n /h 0 ). For T_ < t < T_ + T n 
Brownian motion with a mean square displacement ( x 2 (t )) = 2(ksT/^ n )t causes the liquid-liquid interface to diffuse 
to the next terrace edge, after which thermal motion becomes biased, (x(t)} > 0, and the cycle repeats. 
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CONCLUSIONS 


Theoretical description in the framework of statistical physics predicts that surface nanostructures with directional 
asymmetry can induce nontrivial wetting processes that are beyond the reach of conventional continuum-based models 
(e.g., Lucas-Washburn equation). Fully-atomistic simulations showing close agreement with theoretical predictions 
document the thermally-induced displacement of immiscible fluids confined in a slit channel or pore with a nanoscale 
terraced structure. For the studied nanoscale systems under neutral wetting conditions, a water-air interface would 
exhibit mean displacement rates between 0.1 and 1 m/s. Moreover, the observed fluid displacement can oppose 
the action of wetting or capillary forces. Hence, the studied mechanism can induce the wetting and dewetting of 
nanostructured pores or capillaries under unexpected wettability conditions. The maximum contact angle for which 
thermally driven imbibition or drainage occurs against the action of capillary forces is largely determined by the 
geomtry of the channel and terraced structure. 

The analytical approach presented in this work indicates the necessity to consider the thermal motion of liquid- 
fluid interfaces in order to predict novel phenomena and associated useful effects. The proposed thermal ratchet 
mechanism enabled by engineered surface nanostructures could enable novel nanofluidic devices for passive handling 
and separation. Micro/nanostructured surfaces for superhydrophobic or superoleophobic behavior could exploit the 
thermally driven drainage of micropores and cavities to enhance their performance. In the presence of electrolyte 
solutes, the studied thermal ratchet mechanism could be employed to direct the transport of charges, which can 
potentially enable microfluidic applications for energy harvesting. 
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